load('MG1655_LB_cellmorphology.mat');
temp = MG1655_LB_cellmorphology;
close all;


T = [temp.temperature];
W = [temp.cell_width];
We = [temp.cell_width_err];
L = [temp.cell_length];                 
Le = [temp.cell_length_err];

%% Make violin plot of cell lengths
Lviolin = [];
Tviolin = {};

Ncell = [];
Tcell = {};
for k=1:length(temp)

    lt = temp(k).lengths';
    Lviolin = [Lviolin; lt];
    strtemp = char(strcat(char(num2str(temp(k).temperature)), ' °C'));
    Tcell{k} = strtemp;
    Ncell = [Ncell; length(lt)];

end

T = cell(length(Lviolin),1);
for k=1:length(temp)

    if k==1
        T(1:Ncell(1)) = {Tcell{1}};
    else
        n1 = sum(Ncell(1:k-1));
        n2 = sum(Ncell(1:k));
        T(n1+1:n2) = {Tcell{k}};
    end

end

figure;
violinplot(Lviolin,T)

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
ylabel('Cell lenght (um)')
% xlabel('Temperature (C)')
% ylim([0.8 1.1])
box off;

%% Make violin plots of cell widths
Wviolin = [];
Tviolin = {};

Ncell = [];
Tcell = {};
for k=1:length(temp)

    wt = temp(k).widths';
    Wviolin = [Wviolin; wt];
    strtemp = char(strcat(char(num2str(temp(k).temperature)), ' °C'));
    Tcell{k} = strtemp;
    Ncell = [Ncell; length(wt)];

end

T = cell(length(Wviolin),1);
for k=1:length(temp)

    if k==1
        T(1:Ncell(1)) = {Tcell{1}};
    else
        n1 = sum(Ncell(1:k-1));
        n2 = sum(Ncell(1:k));
        T(n1+1:n2) = {Tcell{k}};
    end

end

figure;
violinplot(Wviolin,T)

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
ylabel('Cell width (um)')
% xlabel('Temperature (°C)')
ylim([0.6 1.3])
box off;

%% Compare cell width distributions using ttest2

pvals_W = zeros(length(W),length(W));

for k=1:length(temp)
    wtest = temp(k).widths;
    for j=1:length(temp)
        wcompare = temp(j).widths;
        [h,p,ci,stats] = ttest2(wtest,wcompare);
        pvals_W(k,j) = p;
    end
end

%% Compare cell length distributions using ttest2

pvals_L = zeros(length(L),length(L));

for k=1:length(temp)
    Ltest = temp(k).lengths;
    for j=1:length(temp)
        Lcompare = temp(j).lengths;
        [h,p,ci,stats] = ttest2(Ltest,Lcompare);
        pvals_L(k,j) = p;
    end
end